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Abstract. In hybrid models, which combine hydrodynamical and transport approaches to describe different 
stages of heavy-ion collisions, conversion of fluid to individual particles, particlization, is a non-trivial 
technical problem. We describe in detail how to find the particlization hypersurface in a 3+1 dimensional 
model, and how to sample the particle distributions evaluated using the Cooper-Frye procedure to create 
an ensemble of particles as an initial state for the transport stage. We also discuss the role and magnitude 
of the negative contributions in the Cooper-Frye procedure. 
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1 Introduction 

In recent years the so-called hybrid models JJH1E1I1MHSH71 
[HHH] have become very popular in describing the expansion 
and evolution of the hot dense matter created in ultrarel- 
ativistic heavy-ion collisions. In these models the different 
stages of the expansion are described using different mod- 
els: The early stages of the expansion, when the matter 
is presumably partonic, are described using ideal or dis- 
sipative fluid dynamics, whereas the late dilute stage af- 
ter hadronization is described using a hadronic transport 
model like RQMD [10], UrQMD [TTj or JAM [T2]. 

Hybrid models are conceptually attractive since they 
attempt to combine the best features of fluid dynamics 
and transport: The complicated microscopic dynamics of 
hadronization can be sidestepped by assuming it to hap- 
pen adiabatically, in which case it can be described using 
fluid dynamics. On the other hand, the dilute hadronic 
matter is assumed to be highly dissipative and deviate 
gradually from local equilibrium. Description of such a 
matter using fluid dynamics is demanding, but is trivial 
using microscopic transport, which by construction can 
handle matter arbitrarily far from equilibrium. Also, the 
transport models describe the last rescattering of parti- 
cles, so-called freeze-out, based on individual scattering 
cross sections. Thus freeze-out is not controlled by a pa- 
rameter which value is known only after comparison with 
data. 

Furthermore, it has been realised that a consistent 
comparison with the data is easier if one generates an 
ensemble of particles event- by- event, and analyses these 
ensembles in a similar way than the actual data has been 
analysed [9l lT3lfT4] . In hybrid models this requires no fur- 
ther effort since the transport models are based on propa- 
gation of individual particles along their scmiclassical tra- 



jectories, and thus require the generation of such ensem- 
bles to begin with. 

However, connecting two approaches with very differ- 
ent degrees of freedom — densities, velocity, and possi- 
ble dissipative currents in the fluid, and individual par- 
ticles in the cascade — is highly nontrivial. Instead of 
running hydro and cascade side by side and using the 
cascade calculation to provide the boundary conditions 
for fluid dynamics, all the present hybrid models solve 
fluid dynamics independently of the cascade taking the 
boundary condition as vacuum at infinity (see discussion 
in Ref. |15j). The boundary where one switches from fluid 
dynamical to transport description is then determined a 
posteriori, once the evolution of the fluid is known. This 
boundary, or switching surface, is usually chosen to be a 
surface of constant temperature, energy density, or time. 
The particle distributions on this surface are evaluated 
using the Cooper-Frye procedure |16| . These distributions 
are sampled to generate an ensemble of particles with well 
defined positions and momenta, and these ensembles are 
used as an initial state for the transport stage. Note that 
even if the transport stage describes the final freeze-out of 
the particles without additional parameters, the choice of 
the criterion where to switch from fluid to cascade is an 
equally free parameter than the freeze-out temperature in 
a conventional hydrodynamical model. Unfortunately the 
final results are also somewhat sensitive to this switching 
criterion as discussed in Ref. [3], and as we will discuss 
later. 

We want to emphasise that this switching from fluid 
dynamics to transport is not freeze-out. By definition, 
there are no rescatterings, only resonance decays after 
freeze-out. Thus, if one switches from fluid to transport 
at freeze-out, the transport stage is unnecessary. It does 
nothing else but lets the resonances decay. On the other 
hand, the switching is not necessarily hadronization either. 
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Since the deconfinement transition is a smooth crossover, 
there is no clear point where hadronization should take 
place, and it is easier to choose to switch at a stage where 
the system has already hadronized. In the following we 
assume that there is no change in the physics nor in the 
properties of the system on the switching surface. It is 
only a change in the description of the system. Thus we 
call the process of changing from fluid dynamics to trans- 
port particlization, conversion of fluid to particles. 

In this paper we describe some technical aspects of 
particlization in hybrid models, and how different choices 
of the switching criterion and constraints imposed to the 
sampling of the distributions affect the final particle distri- 
butions in E lah = 160A GeV Pb+Pb (SPS) and y/s^ = 
200 GeV Au+Au (RHIC) collisions. In particular we dis- 
cuss finding the particlization surface and its properties in 
section [5] the negative contributions of Cooper-Frye pro- 
cedure and how to take them into account when sampling 
the distributions in section [3J the actual sampling pro- 
cedure in section 2] and calculations within one specific 
hybrid approach in section [5] 



2 Surface finding 

The Cooper-Frye procedure [16] for calculating particle 
distributions on a surface is based on evaluating a parti- 
cle four-current through a surface, and a kinetic theory 
decomposition of a four-current j 11 in terms of particle 
distribution f{p,x): 

N = J da^(x) = J da^J ^p»f(x,p) 

= ^da^f(x,p) » £ Ar^/fop). (1) 

Thus one needs to find not only the location of the 
surface a where one applies the Cooper-Frye formula, but 
also its normal. If we knew the analytic expression for the 
surface, its normal would be simply given by |17U18j 

. dx a dx 13 dx 1 . . 

da ll = s ll ^ — — —dadbdc, (2) 

where a, b, c are the coordinates on the surface and £ MCt( 3 7 
is the totally antisymmetric Levi-Civita tensor. However, 
we do not have an analytic expression since we solve the 
evolution of the system numerically. Finding the location 
of an isosurface on a discrete grid is easy, but evaluating 
the siz43 and normal of each discrete surface element is 
not, since one has to make sure that the surface elements 
cover the entire surface without leaving holes or double 
counting any part of the surface. 

In computer graphics and image processing the prob- 
lem of finding an isosurface of a discrete scalar field is 
well known, and many algorithms have been proposed for 
the task, see, e.g. Ref. [TH], and references therein. One 



of the best known of these algorithms is so called March- 
ing Cubes algorithm |20j . a simplified version of which 
was implemented in a hydrodynamical model by Kataja 
and Vcnugopalan [5T], an d subsequently used in AZHY- 
DRO [55]. The original Marching Cubes algorithm was 
extended for hypersurfaces in four dimensional space in 
Ref. |19j . and a simplified version implemented in a 3+1 
dimensional hydrodynamical model by Schcnke |23j . How- 
ever, it is known that the original Marching Cubes algo- 
rithm cannot resolve all possible surface configurations, 
and may leave holes in the final surface [2~ill2"5] . This prob- 
lem is not serious: As quoted in Ref. [53], only 1% of the 
surface elements in a typical heavy-ion collision calcula- 
tion arc not fully resolved. But, if one is doing event-by- 
event hydrodynamical calculations with irregular initial 
conditions, one may expect much more complicated struc- 
tures to appear, and the Marching Cubes algorithm may 
leave more holes in the surface. To avoid this problem alto- 
gether, we have slightly modified the algorithm proposed 
in Ref. |26j . and generalised it for finding a three dimen- 
sional hypersurface in four dimensional spac^H For the 
lack of a better name, we call it here "disordered lines" al- 
gorithm. This algorithm may also have been implemented 
in Ref. [17], but the description there is too vague to tell. 

In algorithms solving the equations of motion of fluid 
dynamics like SHASTA [27], a grid point is considered 
to be in the middle of a corresponding volume element. 
For purposes of surface finding it is useful to consider the 
dual of the grid, where the grid points are thought to be 
at the corners of a volume element, and the values within 
a volume element can be obtained by interpolation. To 
find the location of the surface it is practical to "march" 
through the entire grid and check every volume clement 
whether the values at the corners are all above or all below 
the isovalue. If not, i.e. if some of the corners are above, 
and some below the isovalue, the surface passes through 
this volume element. 

An alternative to this exhaustive search method is the 
continuation method [28) . which after finding a surface 
element, checks only the neighbouring volume elements 
which one of them contains the surface, moves to that 
element, and continues until it either finds the edge of the 
grid or returns to the volume element it found first. The 
continuation method is faster than the exhaustive search, 
but if the surface is disjoint, it will not necessarily find all 
parts of the surface, whereas the exhaustive search method 
by construction finds all parts of the surface. 

In the following we call volume elements (hyper)cubes, 
and explain the algorithm in a 3D case first. Once a cube 
containing the isosurface is found, the positions of the iso- 
surface on the edges of the cube can be found by linear 
interpolation, sec illustration in Fig. [1] How to sort these 
intersection points to form a polygon (or a polyhedron 
in 4D) which does not cross itself is the crucial part of 
any algorithm. It can be done in three different ways |29| : 
By using protomeshes [30], by using a lookup table like 



1 The size enters as the length of the discrete normal vector 
Aa u . 



Fortran and C++ subroutines, Cornelius, implementing 
this algorithm in 3D and 4D, are available at |httpsT7 /karman • 
[physics .purdue . edu/OSCAR/ 



Pasi Huovinen, Hannah Petersen: Particlization in hybrid models 



3 





Fig. 1. A 3D volume element, a cube, with four corners above 
(solid dots) and four below (open dots) the isovalue. In the left 
panel the points where the interpolated value is the isovalue 
are marked with crosses. In the right panel these intersection 
points have been connected to form a polygon. This polygon 
is a surface element of the isosurface. 




Fig. 2. Reduction of a three dimensional problem into a series 
of two dimensional problems: We find the surface element by 
looking for its edges on the faces of the cube. 



in the original Marching Cubes algorithm 19,20,25 , or 
algorithmically [28 l l3T | l32 l [33] like we do here. 

In three or more dimensions the ordering of the inter- 
section points is difficult because of the many possibilities 
involved, but in two dimensions it is almost trivial. There- 
fore we reduce the problem of finding the surface into a 
series of two dimensional problems of finding a line within 
a square, a line which is an edge of the polygon we want 
to find, see Figs. [T]and[2J If all of the corners of a square 
are above or below the isovalue, there is of course no edge 
of the polygon there. If only one of the corners is above or 
below, we have the configuration shown in Fig. [3] a), and 
if two neighbouring corners are above the isovalue, the 
configuration is shown in Fig. [3]b). The case where two 
corners above the isovalue are located diagonally is am- 
biguous. How to choose between configurations shown in 
Figs.[3]c) and d)? Wc follow the original idea of Ref. [26]: If 
the value interpolated at the center of the square is above 
or equal to the isovalue, the center must be inside of the 
surface, and if it is below, the center must be outside of 




Fig. 3. The nontrivial cases of connecting intersection points 
at the edges of the face. 



the surface. A more sophisticated rule has been suggested 
in Rcfs. 25,131]: The values on the face are interpolated 
bilinearly. The isocontours in a bilinear interpolation are 
hyperbolas, and thus one should interpolate the value at 
the center of this hyperbola to decide whether the cor- 
ners above or below the isovalue are linked. A special case 
occurs when the value in any of the cube corners is the iso- 
value. In that case we consider such a corner to be within 
the surface, and place intersection points 10 -9 times the 
edgelength away from the corner, on those edges of the 
cube where the neighbouring corner has a value smaller 
than the isovalue. 

After the intersection points at each face have been 
found, each edge of the polygon is characterised by two 
points: its ends. An important difference between our al- 
gorithm and that of Refs. [21)1151] is that in most common 
cases we do not need to order the edges of the polygon so 
that the next edge in the list begins where the previous 
one ends. We may check the faces of the cube in arbitrary 
order, and keep the edges of the polygon stored in the 
same arbitrary order. Important exceptions to this rule 
are the cases where the face by face search returns six or 
more edges: The surface may consists of two or more dis- 
connected parts, see Fig. @J and we have to sort the edges 
in sequence to find out whether the surface is disjoint, 
and which edges belong to which polygon. We group the 
edges accordingly, and treat the separate polygons inde- 
pendently of each other. 

Note that allowing the surface to consist of several 
parts is necessary for the consistency of the surface. Fig.@] 
depicts when a badly resolved ambiguity on a face 

may lead to a hole on the surface, or counting parts of the 
surface twice. Since we solve the ambiguities by using only 
the information available at the face of a cube, the face is 
always resolved in the same way, no matter in which cube 
it is taken to belong to. Thus the surface elements form 
a consistent surface without holes nor double counted el- 
ements. Furthermore, the rules we described are sufficient 
to resolve any configuration of values at the cube corners. 

After the possible ordering of the edges and dividing 
them in separate groups if they form two or more discon- 
nected surface elements, we evaluate the area and normal 
of the polygon(s). The use of Cooper-Frye procedure re- 
quires not only the surface area and normal, but also the 
values of the hydrodynamical fields on the surface. They 
are best evaluated as interpolated values at the centroid 
of the surface. Since we have to evaluate the centroid for 
this purpose anyway, we triangularise the polygon using 
the centroid: By connecting the ends of each edge of the 
polygon to the centroid, we obtain a set of triangles which 
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double counting 




hole 




Fig. 5. Examples of triangularisation of the polygon in a sim- 
ple and a complicated case. 




consistent 



Fig. 6. A 2D projection of a 4D hypercube. In the figure, 
the cube in the middle is the hyperface of the hypercube at 
coordinate 4=1, the large cube around it is the hyperface 
at t = 0, and the grey lines connecting them are edges with 
constant values of coordinates x,y,z. 



Fig. 4. An example of a configuration where badly resolved 
ambiguity may lead to double counting or a hole on the sur- 
face. The cubes left and right represent neighbouring cubes 
separated for the sake of clarity. 



cover the entire polygon. As depicted in Fig.[5j the polygon 
is not necessarily planar. A sum of the areas and normal 
vectors of the triangles approximates the area and normal 
of the polygon as 



(3) 



where and are vectors from the centroid to the ends 
of the zth edge of the polygon, and /, = ±1 is chosen 
so that each of the normal vectors \fi&i x is directed 
towards lower values, i.e.. outside. 

To make the illustrations of the hypersurfaces in a hy- 
percube more understandable, we first show a simple 2D 
projection of a 4D hypercube in Fig. [51 and how an evolv- 
ing 2D surface in 3D space spans a 3D hypersurface in 
4D spacetime in Fig. [7] This hypersurface element forms 
a polyhedron within a hypercube (Fig. [5]), and like a poly- 



gon can be divided into a group of triangles, a polyhedron 
can be divided into tetrahedra. 

We proceed analogously to the 3D case. We first divide 
the problem into eight three dimensional problems: As a 
cube consists of six squares, a hypercube consists of eight 
cubes, see Fig [9l The surface in each of these is found in 
the same way than described above. The centroid of these 
2D surfaces, which form the faces of the tetrahedra, is 
calculated, and the corners of the triangles forming these 
faces are recorded. This is illustrated in Fig. [5J where the 
triangularisation of the face within the t = hyperface is 
shown. 

As in 3D the ordering of the edges is not important 
unless the surface consists of several disconnected parts. If 
the number of edges and the number of hypercube corners 
above and below isovaluc indicate that this is possible, we 
order the edges to group them according to the polyhedron 
they belong to, and treat the groups as separate surfaces. 

We evaluate an approximative centroid for the polyhe- 
dron. Analogously to the 2D surface in a 3D space, con- 
necting the centroid to the corners of the triangles form- 
ing the faces of the polyhedron, creates a set of tetrahedra 
which fills the volume of the polyhedron, see Fig. [TUJ The 
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t=0 




t=0.2 






t=0.4 



Fig. 8. A hypersurface element, a polyhedron, within a hyper- 
cube. A triangularisation of the polyhedron's face with t = is 
shown. We do not show the triangularisation of the other faces 
for the sake of clarity. 





t=0.6 





t=0.8 





t=l 




Fig. 7. Time evolution of a surface element in four (left) and 
three dimensional (right) representation of a volume element. 
A 2D polygon moving in time spans a 3D polyhedron in a 
hypercube. 



Fig. 9. Reduction of a four dimensional problem into a series 
of three dimensional problems. 



hyperarea, i. e. the volume, of the polyhedron and its nor- 
mal can be calculated as a sum of the volumes and normals 
of its constituent tetrahedra [53] : 



where a*, bi and Cj are vectors from the centroid to the 
corners of the ith triangle, and fi = ±1 chosen so that 
each of the normal vectors of the tetrahedra is directed 
towards lower values, i.e. outside. 

We want to emphasise that this algorithm is by no 
means the only approach for creating consistent surfaces 
with no holes nor double counting [25115X11512] . Especially 
Bernd Schlei's VESTA and STEVE algorithms [30], used 
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Fig. 10. A part of the tetrahedronisation of the polyhedron 
based on triangularisation of the faces of the polyhedron. Full 
tetrahedronisation is not shown for the sake of clarity. 



in Ref. |34j , are faster than ours and produce even higher 
quality surfaces. 

In the recent heavy-ion physics literature two other al- 
gorithms for finding a surface in four dimensions have been 
described. In Ref. [5] the surface elements are taken to be 
the common faces of neighbouring volume elements when 
the isovalue is reached between them. This algorithm is 
simple and robust, but it requires very small grid spacing 
to avoid numerical artefacts (see discussion in Refs. [23J 
[33]). "Very small" in this context means that the volume 
element for the surface extraction should be the same than 
the volume element for the fluid evolution, whereas the 
more sophisticated algorithms discussed here allow signif- 
icantly larger volume elements for surface evaluation than 
for fluid evolution. For example, in Ref. [33] the results 
were calculated using timestep At — 0.04 fm and grid 
spacing Ax = Ay = 0.1 fm, whereas the surface elements 
were evaluated in cubes 10 times larger in time, and 2 to 4 
times larger in both spatial directions (depending on the 
impact parameter). In the following, however, we play it 
safe and use the same grid spacing and timestep for both 
the hydrodynamical evolution and surface extraction to 
calculate the results shown in Sec. [5] 

Another algorithm was described in Ref. [33] and called 
a projection method. To our understanding this method 
insists that there is only one surface element within each 
volume element. As discussed above this may lead to dou- 
ble counting some parts of the surface. 



3 Negative contributions of Cooper-Frye 



One of the conceptual problems in using the Cooper-Frye 
procedure (Eq. ([1])) to describe freeze-out in hydrodynam- 
ical models is that if the surface clement of the freeze-out 



surface is spacclikfl particles with certain momenta count 
as a negative contribution to the distribution: These par- 
ticles move inwards, they are not emitted at the surface 
but absorbccQ. However, in a hybrid model, these nega- 
tive contributions arc only a technical, not a conceptual 
problem. In an ideal case, the switch from fluid to cascade 
is done in a spacetime region where both descriptions give 
(approximately) equal solutions pQ. Thus, in the vicinity 
of the switching surface the particle distributions are close 
to thermal, and correspond to the hydrodynamical solu- 
tion on both sides of the surface. Assuming that the sys- 
tem is dilute enough for the kinetic theory decomposition 
of the fluid variables in terms of particle distributions to be 
valid, the particle distributions f(p,x) are known on the 
surface, and there are particles moving inwards through 
the surface as described by Cooper-Frye. If this is not the 
case, the particlization takes place in a region where trans- 
port and hydro do not give equal solutions, and the switch 
from hydro to cascade should take place elsewhere. 

Thus as an initial state for the cascade, we have to 
count the particles passing through the switching sur- 
face from hydro to cascade. Their distribution is easily 
obtained by augmenting the distribution in Eq. ([!} by a 
(9-function counting particles going outward: 



E 



dN(x) 
dp 3 



da^f(x,p)0(da^ 



(5) 



We have to calculate the number of particles coming from 
hydro to cascade using this distribution, but we must also 
have a drain term in the cascade: All the particles cross- 
ing the particlization surface from cascade to hydro must 
be removed from the cascade. Otherwise the conservation 
laws are not obeyed. It is not possible to remove particles 
from cascade at the rate described by Cooper-Frye pro- 
cedure, but we can calculate the distribution of particles 
which were removed from cascade when they entered the 
fluid dynamical region. This provides a consistency check: 
At the end of the evolution we can compare whether the 
phase-space distribution of particles removed from cas- 
cade is (approximately) equal to the negative contribu- 
tions given by Cooper-Frye. If they are, the calculation 
was consistent, if not, then the switch from hydro to cas- 
cade took place in a wrong place, we have to change the 
switching criterion, and redo the calculation. Of course 
this argument assumes that such a spacetime region ex- 
ists where transport and (dissipative) fluid dynamics give 
equal spacetime evolution. If this is not the case the whole 
raison d'etre of hybrid models is questionable, since the 
switching becomes as arbitrary as the freeze-out in hydro- 
dynamical models. 

Unfortunately the treatment described above is tech- 
nically challenging. First, it is numerically expensive to 



We use the convention to describe a surface as time- or 
spacelike according to its normal vector. dcr M dcr A ' > is time- 
like, and d<7,jdcr M < is spacelike. 

4 In the context of hydrodynamical models, these negative 
contributions have been discussed, and possible cures pro- 
posed, in Refs. [T5ll34U36ll371f38] . 
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evaluate the number of particles emitted when the distri- 
bution is modified by the 0- function (Eq. ([5])). Second, 
we are using UrQMD to describe the hadronic transport, 
and the present version of UrQMD docs not allow tracing 
the propagation of individual particles in such a way that 
whenever a particle crosses the boundary between cascade 
and fluid dynamics, it would be removed from cascade. 
Such a feature would require extensive rewriting of the 
code, and is therefore postponed to later. Furthermore, 
we are using ideal hydro to describe the fluid dynamical 
part, and we do not expect that there exists a region where 
UrQMD would lead to ideal fluid dynamical behavioui0. 

Since the correct treatment is at the time of this writ- 
ing beyond our reach, we have to settle with an approx- 
imative solution like in Ref. [5]: First, we ignore all the 
surface elements where the flow is directed inwards, i.e. 
where der^it^ < 0, since the net number of particles pass- 
ing through such a surface element is negative. There are 
more particles moving inwards than outwards. Second, on 
the surface elements where the flow is directed outwards, 
da^u^ > 0, we use the net number of particles passing the 
surface, dcr M j M , as the number of particles emitted, but we 
use the thermal distribution with (9-function, Eq. as 
their momentum distribution. This leads to a small vio- 
lation of conservation laws. The emitted particles either 
have too large energy, or we generate too few of them, 
depending on the constraints we impose on the sampling 
procedure. 

There arc approaches in the literature where the prob- 
lem of negative contributions is solved differently. Instead 
of sampling the phase-space distributions like we do here, 
in VISHNU [71H0] and in Ref. [13], the momentum distri- 
bution of each particle species is calculated by integrating 
over the particlization surface, Eq. (UJ), like in the con- 
ventional calculation of particle spectra in hydrodynam- 
ical models |411l42] . The number of particles of species i 
per unit rapidity, dNi/dy, in these boost-invariant distri- 
butions is calculated, and the sampling is constrained to 
produce either WdNi/dy particles in a wide rapidity inter- 
val of —5 < y < 5 [4"I] . or dNi/dy particles in an interval 
—0.5 < y < 0.5 [13]. These yields are not integers, how- 
ever. In VISHNU the yield, tOdA^/dy, is rounded down 
to the nearest integer [UJ , whereas in ref. P~3] the decimal 
part is taken as a probability whether to create one addi- 
tional particle over the rounded down yield [12] . At this 
stage the sampled particles carry no information of their 
position. In Ref. |13] there are no rescatterings after par- 
ticlization, and this lack of information is not a problem, 
whereas in VISHNU each particle is given a position ac- 
cording to the probability given by the Cooper- Frye distri- 
bution, neglecting the part where the distribution would 
be negative, Eq. [5] Thus this approach reproduces the 
Cooper-Fryc momentum distribution, but it skews the po- 
sition distribution of the particles reducing the fraction of 
particles emitted from timelike surface elements. How this 
affects the final distribution of particles after rescatterings 
is unknown. Another disadvantage of this approach is that 

5 For the difficulties in approaching the ideal fluid limit using 
transport, see e.g. Ref. [39] . 



since it requires evaluating the momentum distributions of 
all particle species by integrating over the particlization 
surface, it is slow. 

As mentioned in the beginning of this section wc want 
to switch from fluid to cascade in a region where both ap- 
proaches give similar results, and that there is no physical 
change in the region where the switch takes place. This 
means that the solutions to the negative contributions at 
freeze-out proposed in Refs. [Mll5rJll5T] are not applicable 
in our case. In those approaches a non-thermal freeze-out 
distribution which contains no particles moving inwards 
is postulated. But since we assume that the distributions 
are the same on both sides of the switching surface, such 
a distribution on the particle side would be a contradic- 
tion. Furthermore it is perfectly possible that the hadrons 
in the cascade scatter back to the spacetime region de- 
scribed by fluid dynamics, and assuming that there are no 
particles entering that region is an oversimplification. 

4 Sampling the distributions 

Any hadron cascade requires an particle ensemble with 
well defined particle species, momenta, and positions as 
an initial state. To generate such an ensembles we Monte- 
Carlo sample the particle distributions on the particliza- 
tion surface. For this purpose we employ two different 
sampling algorithms. The first one that we call 'allcells' 
sampling, contains a loop over all hypersurface elements 
and in each of the elements, i.e. cells, the particles are 
sampled according to the steps explained below. In this 
case the quantum numbers like energy, net baryon num- 
ber, net strangeness and electric charge are only conserved 
on the average over many sampled events. The other al- 
gorithm that is dubbed 'mode' sampling introduces a way 
to conserve all these quantities event-by-event, as we will 
demonstrate in Section l5T4l The whole sampling algorithm 
is based on previous work published in [5], but now ap- 
plied without the assumption of an isochronous transition 
and with slight improvements. 

In general fluctuations in conserved quantum numbers 
occur in event- by- event studies when only part of the sys- 
tem is described. For example, the spectators are usually 
not included in the hydrodynamical description of the sys- 
tem. As well fluctuations occur in a rapidity interval which 
does not contain all the emitted particles. But, once an 
initial state of the system is defined, and the whole sub- 
sequent evolution is considered, the conserved quantum 
numbers must be conserved. The straightforward 'allcells' 
sampling does not do it, and therefore we have developed 
the more complicated 'mode' sampling approach. 

The first step in sampling hadrons on the hypersurface 
is to decide in which cells a particle has been produced. 
For example, in a central Au+Au collision at the highest 
RHIC energy there are roughly 10 7 hypersurface elements, 
but only ~ 10,000 particles produced. Evaluating first, if 
there is a particle produced in a cell, before doing any of 
the other steps, results in a speed-up of the calculation, 
which is essential for the application to event- by- event cal- 
culations. 
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The number of particles of each hadron species pro- 
duced in one cell is calculated according to the following 
formula (only the surface elements with u^da^ > are 
considered) : 



N, 



(6) 



where rii is the particle density in the local rest frame and 
i is the index of the particle species. It is important to take 
into account all 150 particle species and their antiparticles 
implemented in UrQMD to match the equation of state on 
the particlization surface. Assuming a Boltzmann distri- 
bution the integral over momentum space for the particle 
number density in the local rest frame can be evaluated 
analytically and the result is 



(2tt) 3 



(71 



where <?.; is the degeneracy factor for the respective par- 
ticle species. All the information about the particle prop- 
erties is read in directly from the UrQMD tables to avoid 
any mismatch. The chemical potential includes the baryo- 
chcmical potential and the strangeness chemical potential 
in the following way 



fi = B ■ fi B + S ■ /i S ■ 



(8) 



where S is the quantum number for strangeness and B is 
the baryon number. For pions the Bose distribution has to 
be taken into account because their mass is on the order of 
the temperature of the system. Expanding the distribution 
in series and integrating leads to 



(2tt)* 



fc=i 



kr 



T 



(9) 



where we stop the summation after 10 summands, when 
it has sufficiently converged. 

By summing up all the contributions from different 
species the total number of particles in this cell iV = ^ TVj 
is known. As long as TV < 0.01 one can interpret this 
number directly as a sampling probability and randomly 
decide, if a particle is produced or not. If N is larger, 
one needs to sample a Poisson distribution with N as the 
mean value to decide the actual number of particles pro- 
duced in this surface element. It turns out that the cases 
where 2 or more particles per cell are produced are rare. If 
there is a particle produced the species can be decided by 
sampling according to their probabilities Ni/N. Isospin is 
assigned randomly consistent with the isospin symmetry 
of a system in chemical equilibrium. 

The four momenta of the particles are sampled accord- 
ing to the local Cooper-Frye distribution (only the parts 
where f (x , p)p^ da ^ > are considered) 



(10) 



where f(x,p) are the boosted Fermi or Bose distributions 
corresponding to the respective particle species includ- 
ing again the chemical potentials for baryon number and 



strangeness. Finding the maximum of the distribution and 
then applying the rejection method is crucial since the dis- 
tribution functions in momentum space are highly peaked 
and dependent on the masses and thermodynamic param- 
eters. Currently, an approximative maximum of 



(11) 



is determined by a coarse loop over the three-dimensional 
momentum space. This approximate maximum is multi- 
plied by 1.2 to make sure it is definitely larger than the 
function to be sampled. 

Then, a momentum vector p is chosen randomly, and 
an additional random number between zero and the above 
mentioned maximum of the distribution is generated. If 
this random number is smaller than the value of the dis- 
tribution at this momentum, the momentum is accepted. 
If not, another momentum and another random number 
are generated, and the process is repeated until an accept- 
able momentum is found. 

As described in Sec. [3] we neglect negative parts of 
the distribution functions which slightly alters the result- 
ing rapidity and transverse momentum spectra. It is im- 
portant to sample the momenta according to the boosted 
distribution function instead of sampling the equilibrium 
distribution in the local rest frame and then boost the 
momentum four-vector back to the computational frame. 
The second procedure leads to a violation of energy con- 
servation since one does not reproduce the whole tensor 
T^ v in the computational frame correctly |16j . 

Imposing strict conservation laws on grand canonical 
distributions on event-by-event basis is not wholly consis- 
tent, and we minimise any bias it creates by the following 
procedure: To conserve energy, net baryon number, net 
strangeness and electric charge in all events we do seven 
subsequent random loops called modes over the hyper- 
surface. During the first mode, cells are randomly chosen 
and particles are sampled until the total energy is con- 
served. From this set of particles only the ones containing 
a s quark are kept, and the rest are discarded. In the 
second mode we produce particles in a similar way, until 
we've produced as much anti-strangeness as the first mode 
produced strangeness. We keep the anti-strange particles, 
and discard the rest. In modes three and four we repeat 
the same procedure keeping the non-strange baryons and 
non-strange antibaryons, respectively, but requiring that 
the net baryon number of these particles is the net baryon 
number of the system. Modes five and six take care of 
the conservation of electric charge by keeping the neg- 
ative and positive non-strange mesons, respectively, and 
finally mode seven takes care of the energy conservation 
by sampling neutral non-strange mesons until the energy 
of the particles corresponds to the energy of the fluid. 

5 Results 

5.1 Description of the Hybrid Model 

The hypersurface finding and sampling algorithms described 
in Sections [2] and U have been implemented in the hybrid 
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model of Rcf. [5] . We show here some preliminary results 
to demonstrate some general trends and overall behaviour 
of the improved model. Note that we have not tried to 
fine-tune any parameters to describe experimental data 
at this point. 

Smooth initial conditions for the hydrodynamic evolu- 
tion have been generated by averaging over 100 UrQMD 
events that are run up to the starting time of i s tart = 2.83 
fm for Pb+Pb collisions at E uh = 160A GeV (SPS) 03] 
and istart = 0.5 fm for Au+Au collisions at £ cm = 200 A 
GeV [II] (RHIC). 

On the constant time surface t = t star t the energy, mo- 
mentum and net baryon number densities are determined 
by representing each particle in the UrQMD event with a 
three-dimensional, longitudinally Lorentz-contracted, Gaus- 
sian density distribution of thermalised matter, which has 
the same energy, momentum and baryon number than the 
particle it represents. Summing up the distributions rep- 
resenting the particles, we obtain a distribution of matter 
in thermal equilibrium, which has the same energy, mo- 
mentum and baryon number than the UrQMD event in 
question. For the RHIC initial state we consider only par- 
ticles within the rapidity interval — 2 < y < 2, whereas 
we include all particles that have interacted at least once 
when evaluating the initial state for SPS. The mean value 
of these distributions from 100 events leads to a smooth 
profile but still includes finite initial velocities in all three 
directions [15] . 

The (3+l)D ideal hydrodynamic evolution is solved in 
Cartesian coordinates using the SHASTA algorithm 05] 
14"?] . The equation of state is calculated within a chiral 
model coupled to the Polyakov loop to include the decon- 
finement transition that reproduces nuclear ground state 
properties and lattice QCD data at zero net-baryon chemi- 
cal potential [48] . This equation of state is also applicable 
at finite net baryon densities, which allows us to calcu- 
late heavy ion collisions at SPS energies to explore the 
beam energy dependence of our findings. In the hadronic 
phase this equation of state is equivalent to the effective 
equation of state of UrQMD, and has the same degrees 
of freedom than UrQMD. The conservation laws are thus 
automatically fulfilled on particlization surface when we 
use Cooper-Frye prescription to calculate the particle dis- 
tributions. 

The hypersurface finder 'cornelius' has been imple- 
mented in the hydrodynamic code in such a way that 
the hypersurface for any transition criterion is evaluated 
while the hydrodynamic evolution is calculated. Then we 
apply the two sampling algorithms and further calculate 
resonance decays using UrQMD. To obtain final results 
that can be compared to experimental data the hadronic 
transport approach is used in addition to calculate the 
rescatterings. 

5.2 Structure of the Hypersurface 

The equations of motion of relativistic hydrodynamics are 
nothing more than an application of conservation laws for 
energy, momentum and charge(s). Thus one of the first 



Table 1. The conservation of energy (E) and net baryon num- 
ber (B) in Au+Au collisions at E cm = 2004 GeV. The values 
in the final state are split into two parts: pos. is flow through 
elements where the energy or baryon flow is directed outwards, 
dcr M T A ' > or dcr M n_Bit M > 0, respectively, whereas neg. is flow 
through elements where energy or baryon flow is directed in- 
wards, da^T^ < or da^UBU 11 < 0, respectively, and thus 
counts as negative. Note that these are not the negative con- 
tributions of Cooper-Frye, see the text. The upper two rows are 
for central (6 < 3.4 fm) and the lower two rows for mid-central 
(6 = 7 fm) collisions. 





E [GeV 




B 




total 


pos. 


neg. 


total 


pos. 


neg. 


initial 
final 


5431 
5430 


5861 


-431 


93.23 
92.74 


97.74 


-5.00 


initial 
final 


2327 
2336 


2455 


-119 


35.84 
35.80 


37.10 


-1.30 



Table 2. The conservation of energy (E) and net baryon num- 
ber (B) in Pb+Pb collisions at E^b = 160A GeV. The layout 
is the same as in Table [T] 





E [GeV 




B 




total 


pos. 


neg. 


total 


pos. 


neg. 


initial 
final 


3117 
3120 


3208 


-88.5 


347.6 
345.9 


355.8 


-9.9 


initial 
final 


1528 
1528 


1570 


-41.5 


170.23 
169.22 


174.26 


-5.04 



checks on the accuracy of the numerical solutions of hy- 
drodynamics is to confirm the validity of the conservation 
laws. 

We show the total energy (E) and net baryon number 
(B) in collisions at two different centralities at RHIC and 
SPS in tables [T] and [2j respectively. The initial value is 
evaluated in the beginning of the hydrodynamical evolu- 
tion by summing over all fluid cells within the particliza- 
tion surface. With the switching criterion we use here, 
e = 2eo, where eo refers to the nuclear ground state en- 
ergy density of 146 MeV/fm 3 , this means summing over 
the cells with e > 2eo- Since the final state we are inter- 
ested in is the isosurface with e = 2eo, we evaluate the 
final state energy and net baryon number as energy and 
baryon number flows through this surface: 

E = T^ddf, and B = nsu^da^, (12) 

J a Jo 

where the energy momentum tensor has been evaluated 
as T» v = (e + P)u^u y - g^P. As seen in the tables both 
the energy and baryon number are nicely conserved with 
better than 0.6% accuracy in all cases. 

However, a closer look at the properties of the parti- 
clization surface reveals that there are surface elements 
where the flow is directed inwards, da fJ ,u tl < 0. As de- 
scribed in Sec. 2] such elements cannot be used as a source 
for sampled particles, and thus we check how large uncer- 
tainty these elements cause to the total energy and baryon 
number. In tables [TJ and [21 the flows through the surface 
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Fig. 11. The position of the e = 2eo isosurface in zt-plane at 
x — y = 0. The left side of the plot (SPS) depicts the surface 
in central Pb+Pb collisions at E\ a \y = 160A GeV and the right 
side (RHIC) in central Au+Au collisions at E cm = 200,4 GeV. 
The dashed lines depict regions where the flow may be directed 
inwards on some surface elements. 



are also separated in two parts: pos. depicts the energy and 
baryon flows through surface elements where da^T^ > 
for energy and da^nBU^ > for baryon flow, whereas 
neg. is the flow through elements where da^T^ < for 
energy and dcr^ngu^ < for baryons. As seen the surface 
elements where flow is directed inwards cause a 3-8% un- 
certainty in the total energy, and a 3-5% uncertainty in 
the total baryon number of the system. Note that because 
of the pressure term in T 00 it is possible that we have an 
element where energy flow is directed inwards but baryon 
flow outwards, or vice versa. 

To understand this phenomenon we located the surface 
elements with da^u^ < and found out that they are al- 
most entirely located at the edges of the system where the 
fireball expands longitudinally, see Fig. [TTJ and there are 
none at midrapidity. Furthermore, the regions depicted by 
dashed lines in Fig. [TT] are not regions where the flow is 
directed inwards, but they are regions where surface ele- 
ments with inwards directed flow randomly appear among 
elements where the flow is directed outwards. These are 
the regions where the longitudinal pressure gradient is still 
large, the fluid flow is very large, v z ~ 0.9-0.95, and the 
isosurface moves outwards with a comparable speed as 
well. It is known that this is a difficult terrain for SHASTA 
and for any algorithm to solve accurately |46j . Since the 
flow velocity is almost aligned with the surface, it does not 
require a large numerical error to flip an outwards flow to 
inwards. Thus it is possible that these elements with flow 
directed inwards are just numerical artefacts. However, it 
is worth remembering that in event-by-event calculations 
the initial states are highly irregular, and occasionally flow 
can be directed inwards for physical reasons. 

Since the common feature of the elements with inwards 
directed flow is the large longitudinal flow velocity, we 
further characterise the surface by it. We bin the surface 
elements according to the magnitude of the longitudinal 
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Fig. 12. The magnitude of the fluid flow through different 
parts of the particlization surface at different longitudinal flow 
velocities in central Pb+Pb collisions at £i a b = 160A GeV 
(SPS) and in central Au+Au collisions at E cm = 200A GeV 
(RHIC), see the text. 



flow velocity \v z \, and integrate over the elements in each 
bin to obtain the magnitude of fluid flow through sur- 
face, |d<T M u M |, at different velocities \v z \. The results are 
shown in Fig. 1121 The contribution from elements with 
flow directed inwards is depicted by the thick curve, and 
as claimed it is concentrated at large velocities. Further- 
more it can be seen that the fluid flow inwards is small 
compared to the fluid flow outwards, which is depicted by 
the thin solid and dashed histograms in Fig. [T2] The 5- 
8% uncertainties in energy and baryon number listed in 
tables [T] and [2] arc thus concentrated at large longitudi- 
nal flow velocities, and can be expected to have only a 
tiny contribution to observables at midrapidity. The rea- 
son why there is a difference between fluid flow da^u 11 
and energy and baryon number flows, da^T^ and da^ j^, 
respectively, is that the former has an additional gamma 
factor and a pressure term, whereas the latter is also sen- 
sitive to the baryon chemical potential fis which is not 
uniform on the surface. 



5.3 Negative contributions 

It is worth remembering that the negative contributions 
to energy and baryon number discussed above in Sec. 15.21 
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Fig. 13. The ratio of the negative contribution to the total 
thermal pion rapidity distribution in central Pb+Pb collisions 
at _Ei a b = 160^1 GeV (SPS) and in central Au+Au collisions at 
E cm = 200,4 GeV (RHIC). The curves correspond to different 
kinds of contributions, see the text. 



are not the same thing than the negative contributions of 
Cooper- Frye discussed in the literature (see Sec. |3j). The 
former are caused by the collective flow being directed in- 
wards, whereas the negative contributions of Cooper-Fryc 
are caused by individual particles moving inwards through 
the particlization surface. As discussed in Sec. SI the sam- 
pling routine cannot take into account particles moving 
inwards nor anything emitted on surfaces where flow is 
directed inwards. Therefore we evaluated these negative 
contributions to see how much omitting them distorts the 
spectra. 

We evaluated the spectra of thermal pions, kaons and 
protons on the particlization surface using the conven- 
tional Coopcr-Fryc procedure, Eq. (|T|), i.e. we integrated 
over the surface to obtain distributions, but we divided 
the contribution into four parts according to whether the 
flow was directed in- or outwards (da^u^ < or > 0) 
and whether the particles with the momentum in ques- 
tion were heading in- or outwards (dcr^p^ < or > 0). It 
turned out that in general the negative contributions to 
pion distribution are largest, so we concentrate on them. 
Some of the actual spectra are shown in Sec. 15.51 where 
we compare the sampled distributions to the integrated 
ones. Here we emphasise the size of the negative contri- 
butions by showing their relative contribution to the total 
spectrum. 



In Fig. 1131 different negative contributions to the ther- 
mal pion rapidity spectrum are shown for collisions at 
RHIC and SPS. To check how much uncertainty the con- 
tribution from elements with flow directed inwards cause, 
we have calculated the contribution from them (dashed 
line, da^u^ > 0) even if contribution from them is not 
negative in the Coopcr-Frye sense. It is seen that at midra- 
pidity contribution from them is 2-5%. A modest contri- 
bution at midrapidity could be expected, since as shown 
in Fig. [T2l the longitudinal flow velocity on these elements 
is large. Interestingly the contribution from these elements 
at large rapidities is very different at RHIC and SPS: The 
magnitude of the relative contribution decreases with in- 
creasing rapidity at SPS, but increases at RHIC. Again 
a result of different flow and emission patterns shown in 

Fig.m 

We may take the surface and flow pattern as granted, 
and simply evaluate the fraction of inward moving par- 
ticles. This is depicted by the solid line, da^p 11 < 0, in 
Fig. 1131 At midrapidity the contribution of such a parti- 
cles is surprisingly large, 10-15% with a larger contribu- 
tion at SPS than at RHIC. That the negative contribu- 
tion is larger at SPS was already implied in Fig. [T2"l where 
the fluid flow through the particlization surface was di- 
vided into a flow through time- and spacelike surfaces. At 
SPS the relative fraction of the fluid flow passing through 
spacclikc surfaces is larger than at RHIC. Since negative 
contributions appear only on spacelike surfaces, the flow 
pattern at SPS provides a possibility for larger negative 
contributions than at RHIC. Furthermore, if surfaces are 
similar, the lower the flow velocity, the larger the negative 
contributions. At SPS the flow velocity is lower than at 
RHIC, and thus it is rather surprising that the difference 
in negative contributions is not larger than what is shown 
in Fig. [lj 

Finally, we want to check how large contribution our 
sampling routine misses. Since it cannot sample the par- 
ticles going inwards nor particles emitted when the flow 
is directed inwards, sampling misses both contributions. 
This is shown as the dotted line in Fig. [13] At midrapidity 
the missed part is totally dominated by the conventional 
negative Cooper-Frye contribution. Thus it is safe to say 
that what comes to sampling particles at midrapidity, the 
inwards directed flow is not a problem, but the conven- 
tional negative contributions of Cooper-Frye can be. 

To further study the effect of negative contributions, 
we show the contributions to px distribution of thermal 
pions at midrapidity, — 1 < y < 1, in Fig. 1141 Again we 
see that the contribution from the inwards directed flow 
is smaller than from inwards moving particles. The lat- 
ter can even reach 40% at low values of transverse mo- 
mentum! Fortunately the size of the negative contribution 
decreases rapidly with increasing transverse momentum, 
and already around 500 MeV they are much more accept- 
able 8-10%. That the negative contributions are largest at 
small px is understandable since the high px particles are 
mostly emitted in regions where the flow velocity is large 
and parallel to momentum, whereas the negative contri- 
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Fig. 14. The ratio of the negative contribution to the total 
thermal pion p-r-distribution at midrapidity —1 < y < 1 in 
central Pb+Pb collisions at £ lab = 160,4 GeV (SPS) and in 
central Au+Au collisions at E cm = 200,4 GeV (RHIC). The 
curves correspond to different kinds of contributions, see the 
text. 



butions arise when the momentum is antiparallel to flow 
velocity. 

We also evaluated the pr averaged «2 of thermal pions 
at particlization in mid-central (6 = 7 fm) collisions at 
RHIC and SPS. When t>2 was evaluated in a conventional 
fashion, we got v 2 = 0.064 and 0.066 at RHIC and SPS, 
respectively. When we removed the negative contributions 
from the distribution, w 2 at RHIC and SPS decreased to 
0.058 and 0.06, respectively. If we also removed the re- 
maining positive contribution from elements with flow di- 
rected inwards, further change on i>2 was on the level of 
10~ 4 . This additional uncertainty of 10% should be kept 
in mind when discussing the elusiveness of the QGP shear 
viscosity |49j . 



5.4 Conservation of Quantum Numbers 

In Figs. [15] and [TB] we show the probability distributions 
of relevant quantum numbers using the two sampling al- 
gorithms described in Section @] As explained, the 'mode' 
sampling is constructed to obey conservation laws event- 
by-event, whereas in the simpler sampling of all hypersur- 
face elements the conserved quantities fluctuate around 
the values given by the positive contributions to the spec- 
tra. Note that these values are not those given in Tables [TJ 
and[2]as energy and baryon number flow through the parts 
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Fig. 15. Probability distribution for the total energy of an 
event in an ensemble of 500 events created using two different 
sampling algorithms in central Au+Au collisions at E cm = 
200/1 GeV. 



of surface where flow is directed outwards, but by the dis- 
tributions of particles coming out through these parts of 
the surface, Eq. (|5|), after these distributions have been 
scaled to yield the same number of particles than the local 
net flow of particles through the surface element, da^nu 11 . 
Therefore, we expect to see a higher mean value in the 
energy and net baryon number in the allcells sampling 
compared to the mode setup. One of the reasons to im- 
plement the mode sampling was to take into account the 
distortions caused by different kinds of negative contribu- 
tions by enforcing the conservation lawtQ. 

The probability distributions for the different quanti- 
ties are rather wide in the allcells sampling, so individual 
events in the more often applied algorithm can have quan- 
tum numbers that are far away from their actual values 
in that event. In a sense this is reasonable since the equa- 
tion of state used during the hydrodynamical stage, and 
the particle distributions sampled at particlization assume 
a grand-canonical ensemble where conservation laws are 
obeyed only in average. But if we aim at a description of 
the collisions on an event-by-event basis, we have a contra- 
diction: In nature conservation laws are obeyed in every 
single event. Thus it makes sense to require that the con- 
servation laws are strictly obeyed during all stages of the 
model, and check how much the fluctuations created by 
the sampling procedure affect the observables. 



6 Remember that in event-by-event simulations inwards flow 
can be physical. 
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Fig. 16. Probability distribution for the total net baryon num- 
ber (B), net strangeness (S) and net charge (CH) of an event in 
an ensemble of 500 events created using two different sampling 
algorithms in central Au-f Au collisions at E cm = 200^4 GeV. 



5.5 Tests of Sampling Algorithm 

First of all, we need to check, whether the sampling algo- 
rithm reproduces the numbers of particles for each species. 
In Figs. [T7] and [T5] the multiplicities of selected particle 
species are compared to the integrated results. Integra- 
tion means here summing up as defined in equation [6] 
for all hypersurface elements where da^u^ > 0. All the 
results in this section are comparisons of thermal yields 
for individual species only, the resonance contribution has 
not been included. 

We have plotted in Figs. [T7] and [TS] the relative de- 
viations from the integrated result, sampled-^ inte- 
grated)/^ integrated, for the two sampling algorithms. 
The deviations from zero are within the expected statis- 
tical fluctuations for 500 events taking into account their 
individual abundances for the allcells sampling. Overall, in 
the allcells sampling the particle multiplicities are nicely 
reproduced at RHIC and SPS energies, whereas the mode 
sampling leads to slightly less particles at RHIC. This re- 
sult is expected, since we enforce a lower value for the 
total energy of the system. The fluctuations of the par- 
ticle yields around that lower value are again purely sta- 
tistical. Since at SPS the inwards directed flow is smaller 
this difference is also smaller and therefore the difference 
between allcells and mode sampled yields is smaller. As 
one can see in the rapidity spectra (Figs. [P31 and EH"]) , at 
SPS energies mode sampling produces more particles than 
allcells sampling. 

The next step is to compare the distribution of the par- 
ticles in momentum space in terms of rapidity and trans- 



Fig. 17. Relative difference of sampled and integrated particle 
multiplicities in central Au+Au collisions at E cm = 200^4 GeV. 
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Fig. 18. Relative difference of sampled and integrated particle 
multiplicities in central Pb+Pb collisions at Ei^b = 160^4 GeV. 



verse momentum. The different lines in Fig. IT9l correspond 
to the integrated result, which is split up into the differ- 
ent contributions depending on the direction of flow and 
momentum with respect to the particlization surface as 
described in the legend, sec Sec. 15.31 

For clarity all the contributions are displayed only for 
pions while for kaons and protons only the full result and 
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Fig. 20. Comparison of the integrated and sampled trans- 
verse momentum spectra for pions, kaons and protons in cen- 
tral Au+Au collisions at E cm = 200A GeV. 



Fig. 19. Comparison of the integrated and sampled rapidity 
spectra for pions, kaons and protons in central Au+Au colli- 
sions at E crn = 200A GeV. 



the result from hypcrsurface elements with flow directed 
outwards are shown. The full line corresponds to the full 
result from the complete hypersurface integration that one 
ideally wants to reproduce. The dashed line shows the re- 
sult for the contribution from the cells where the flow is 
directed outward; these are the cells we take into account 
when sampling the particle number densities. The momen- 
tum sampling on the other hand disregards all the negative 
parts of the distribution function. This treatment slightly 
skews the spectrum such that we end up with a result 
that is in between the dashed and the dot-dashed line for 
the sampled spectra. It may look surprising that the sam- 
pling does not reproduce the yield from outward flowing 
elements at midrapidity. After all, the sampling procedure 
has been normalised to reproduce the yield from those el- 
ements. The reason for this deviation is that the system in 
our calculations is not boost invariant. If it were, the edges 
of the system in longitudinal direction would be infinitely 
far away, and they would not contribute to midrapidity. 
But in our non-boost invariant system they are close, and 
we see the same to happen in longitudinal direction than 
what is depicted in Fig. [14] for transverse direction: The 
low momentum, i.e. small rapidity, region is depleted be- 
cause of the negative contributions. The distributions we 
sample do not contain this depletion, and thus we end up 
having too many particles at midrapidity. On the other 
hand, since the sampling procedure reproduces the total 
yield of particles as shown in Figs. [171 andlTSl we must have 
fewer particles somewhere. A closer look at the spectra re- 
veals that the sampled yield is below the der^u^ > curve 
at rapidities 1 < \y\ < 2.5, and this deficit is sufficient to 
cover the excess at midrapidity. 
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Fig. 21. Comparison of the integrated and sampled rapidity 
spectra for pions, kaons and protons in central Pb+Pb colli- 
sions at _Bi a b = 160A GeV. 



The transverse momentum spectrum is shown in Fig. 
I2TJ1 Due to the logarithmic scale in this case, the differ- 
ences shown in Fig [T?] disappear, and only a small dif- 
ference between the elements with flow directed outwards 
and the full result is visible at low transverse momentum. 
The sampled results are in very good agreement with the 
integrated results. All possible deviations (apart from sta- 
tistical noise) are smaller than the size of the symbols. 

The comparison of rapidity spectra and transverse mo- 
mentum spectra at SPS, as shown in Figs. P2"T1 and P2"2l leads 




to very similar results. The only difference is that in the 
mode sampling there are more pions produced than in the 
allcclls sampling as mentioned above. The reason for this 
has to be investigated further in the future. 

Overall, the sampling algorithms reproduce the parti- 
cle yields and spectra very well within the expected devia- 
tions from the assumptions that are made in the different 
algorithms. The negative contributions are also smaller 
than 5 % in the regions where most of the particles are 
emitted and affect heavier particles less than lighter ones. 



Fig. 23. Yields at midrapidity (\y\ < 0.5) as a function of 
switching criterion in multiples of eo for four different particle 
species (7r~, K + , P, A) in central (6 < 3.4 fm) Au+Au collisions 
at ^/sNN = 200 GeV. The symbols indicate the result with 
resonance decays ('Hydro') and the lines show the result after 
rescattering ('After'). The black stars show the result from 
the previously used 'gradual' transition. The bands on the left 
indicate data by the STAR and BRAHMS collaborations [53 



5.6 Effects of Rescattering 

In this Section we present selected results for multiplic- 
ities, spectra and elliptic flow to demonstrate the effect 
of the hadronic rescattering and the switching criterion. 
All the results are calculated using the hybrid approach 
described in Sec. I5.1l with mode sampling enforcing exact 
conservation laws. To obtain the results directly at the 
switching hypersurface we run UrQMD for 0.1 fm and use 
it only to calculate the resonance decays ('Hydro'). The 
second option is to run UrQMD for 200 fm and include 
all the rescattering dynamics in the hadron transport ap- 
proach ('After'). 

In Fig.[22]thc yields of pions, kaons, protons and Lamb- 
das are shown in central heavy ion collisions at the highest 
RHIC energy The final results using allcclls sampling are 
within 3% to the ones shown here (see Fig.[T7]). Therefore, 
we decided to restrict the number of shown results to mode 
sampling only. The switching criteria that we have chosen 
correspond to the following temperatures 3eo « 163 MeV, 
2e w 154 MeV and 1.5e « 149 MeV and are in the ball- 
park of switching criteria that are commonly used in other 
hybrid approaches [Tl|2l|3l|il|5T[6l[7Pl|g] . The energy density 
eo refers to the nuclear ground state energy density 146 
MeV/fm 3 . The switching criterion needs to be adjusted to- 
gether with equation of state and parameters defining the 



initial conditions to achieve a good agreement with experi- 
mental data. In this analysis we concentrate on presenting 
general features of the results and have not searched for 
the best agreement with the measurements yet. 

The yields are slightly higher for higher switching en- 
ergy densities, because in equilibrium at this temperature 
range, higher temperatures lead to higher yields. The devi- 
ations from equilibrium caused by the sampling procedure 
are not large enough to break this rule. The rescattering 
leads to a reduction of the kaon and proton yield, whereas 
pions and Lambdas are not so much affected. The kaon 
and proton yields are also higher in the iso-energy density 
transition scenario than in the previously employed grad- 
ual transition, where full transverse slices are sampled on 
an isochronous surface, when the whole slice has diluted 
below 5eo- Apart from the fact that the kaon yield is higher 
than the experimental data our results are in reasonable 
agreement with the data. 

The mean transverse momentum as a measure of the 
radial flow development is shown in Fig. [Ml In this case 
the dependence on the switching criterion is opposite to 
the one of the yields. There are two reasons for this be- 
haviour: The first is that for lower switching energy densi- 
ties the system stays longer in the hydrodynamic evolution 
and the particles gain more transverse flow. The second 
argument is that the total energy needs to be conserved 



16 



Pasi Huovinen, Hannah Petersen: Particlization in hybrid models 



1.4 



1.2 



> 

CD 

2. 1.0 

E 
o 

A 
CL 

V 



• if, Hydro 

■ K + 

♦ P 

A A 

After 



h 1 1 1 1 1 r 




0.8 



0.6 



0.4 



▲ 
♦ 



-J i I i I i I i I i l_ 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.! 

X*e 



Fig. 24. Mean transverse momentum at midrapidity (\y\ < 
0.5) as a function of switching criterion in multiples of eo for 
four different particle species (ir~, K + , P, A) in central (b < 3.4 
fm) Au+Au collisions at -^/snn = 200 GeV. The symbols indi- 
cate the result with resonance decays ('Hydro') and the lines 
show the result after rescattering ('After'). The black stars 
show the result from the previously used 'gradual' transition. 
The bands on the left indicate data by the PHENIX collabo- 
ration [54]. 



and if there are more particles produced, there is less en- 
ergy remaining to give them kinetic energy in terms of 
transverse flow at the switching transition. Depending on 
their hadronic cross-sections and their mass the particles 
get pushed a lot (e.g. protons) or almost not at all (pions) 
during the hadronic rescattering stage. 

The transverse mass spectra at RHIC and SPS (shown 
in Figs. [25] and [2S1 respectively) confirm this behaviour. 
The pion spectra are almost identical before and after the 
hadronic cascade, but the kaon and proton spectra begin 
to resemble the experimental data only after the rescat- 
tering. The comparison to the previously imposed gradual 
transition shows that a true iso-energy density switching 
criterion improves the slope of the spectra at high trans- 
verse masses drastically. This can be easily understood, 
since in the gradual transition scenario the full slice needs 
to reach the energy density criterion which delays the 
switching to hadron cascade for the edges in the transverse 
direction. The edges gain very large transverse flow veloc- 
ity, which makes the distributions flatter, and even if the 
edges get very cold, the yield of heavy particles does not 
drop so much that the effect of large flow velocity would be 
negated. What comes to the pT-chstributions, the hybrid 
model based on ideal hydrodynamics works better at high 
RHIC energies than at high SPS energies. This indicates 
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Fig. 25. Transverse mass spectra at midrapidity (\y\ < 0.5) 
for pions, kaons and protons in central (b < 3.4 fm) Au+Au 
collisions at ^/snn = 200 GeV. The dashed lines indicate the 
result with resonance decays ('Hydro') and the full lines show 
the result after rescattering ('After'). The black dotted line 
represents the result from the previously used 'gradual' tran- 
sition. The symbols indicate experimental data by the STAR, 
PHENIX and BRAHMS collaborations [55U53U54TT56] . The pro- 
ton data by STAR has been multiplied by 0.6 to correct for the 
feed-down from Lambdas. 
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Fig. 26. Transverse mass spectra at midrapidity (\y\ < 0.5) 
for pions, kaons and protons in central (b < 3.4 fm) Pb+Pb 
collisions at i5i a b = 160^1 GeV. The dashed lines indicate the 
result with resonance decays ('Hydro') and the full lines show 
the result after rescattering ('After'). The black dotted line 
represents the result from the previously used 'gradual' tran- 
sition. The symbols indicate experimental data by the NA49 
collaboration I57II58I. 
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Fig. 27. Elliptic flow of pions as a function of transverse mo- 
mentum at midrapidity (\y\ < 0.5) in mid-central (6 = 7 fm) 
Pb+Pb collisions at Ei^b = 160^4 GeV. The dashed lines in- 
dicate the result with resonance decays ('Hydro') and the full 
lines show the result after rescattering ('After'). The black dot- 
ted line represents the result from the previously used 'grad- 
ual' transition. The symbols indicate experimental data by the 
NA49 collaboration [59] . 



Fig. 28. Elliptic flow of protons as a function of transverse 
momentum at midrapidity (\y\ < 0.5) in mid-central (b — 7 
fm) Pb+Pb collisions at E lah = l&GA GeV. The dashed lines 
indicate the result with resonance decays ('Hydro') and the 
full lines show the result after rescattering ('After'). The black 
dotted line represents the result from the previously used 'grad- 
ual' transition. The symbols indicate experimental data by the 
NA49 collaboration [59] . 



that the non-equilibrium dynamics gains importance at 
lower beam energies. 

In Figs. Wf\ and P251 elliptic flow as a function of trans- 
verse momentum is shown for pions and protons at high 
SPS energy. Quite surprisingly the pr-differential elliptic 
flow is quite insensitive to the switching criterion. It is 
known that in ideal fluid hydrodynamics pion i>2(pt) ls 
quite insensitive to the freeze-out temperature [35] . We 
presume that the same holds also when one uses the grad- 
ual criterion instead of constant temperature/density for 
particlization, and thus the elliptic flow of pions at par- 
ticlization is very similar in both cases, and evolves simi- 
larly during the cascade. Thus the final V2{px) is similar as 
well. The proton V2(px) shows more sensitivity than pion 
v 2(pt), but the difference is at low px, not at high px 
which is influenced by the different evolution of the edges 
as explained when discussing the px spectra. On the other 
hand, the lower proton vi at low px is in line with the 
ideal fluid expectations where lower freeze-out tempera- 
ture leads to lower proton v^ijpx) at low px- Thus the low 
px range of protons is mostly influenced by the center of 
the system which evolves hydrodynamically much longer 
when one uses isodensity switching than when one uses 
gradual switching. When the hadronic cascade rescatter- 
ing is taken into account, the proton flow at low transverse 
momenta seems to even turn negative as recently observed 
by the CERES collaboration [HQ]. 

6 Conclusions 

In this paper we have discussed a crucial part of hybrid 
models in detail: How to switch from hydrodynamical de- 



scription to cascade. We have described an algorithm to 
find an isovalue surface where the transition takes place, 
the unavoidable negative contributions of Cooper- Frye pro- 
cedure on such a surface, and how to make the sampling 
algorithm, which creates the initial state for the cascade, 
to work on an arbitrary surface and deal with the negative 
contributions. We have seen that in realistic calculations 
the negative contributions are not a large problem, but 
they create an uncertainty of their own which should be 
kept in mind when one draws conclusions from the results 
of the present hybrid or hydrodynamical models. In the 
long run we hope to develop a model where these negative 
contributions are properly treated, and not just ignored. 

We think that the machinery we have described is par- 
ticularly suitable for a proper event- by- event analysis of 
heavy-ion collisions. In event-by-event calculations the ini- 
tial state of hydrodynamics varies wildly from event to 
event, and thus one may expect the particlization sur- 
face to have a very complicated structure. This poses no 
problem for the algorithm described here since it can cre- 
ate consistent surfaces without holes nor double count- 
ing for any configuration. One of the reasons for study- 
ing heavy-ion collisions event-by-event are fluctuations. 
In those studies it is important to distinguish at which 
stage of the evolution fluctuations are created and how. 
The sampling algorithm we described provides an impor- 
tant check by allowing strict conservation of energy and 
charges thus guaranteeing that their fluctuations are not 
caused by the sampling procedure. 

All this machinery is not limited to switching from 
ideal fluid to cascade, but it can be applied at the inter- 
face of viscous hydrodynamics and cascade as well. All 
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one needs to do is to modify the particle distributions at 
particlization accordingly. 

However, we presented here preliminary results only 
without attempting to describe the data well. The next 
step is to identify the main parameters of the hybrid calcu- 
lation and perform a multi-parameter analysis compared 
to bulk observables at RHIC using a sophisticated statis- 
tical emulator to limit CPU time. Once a good switching 
criterion has been identified, the beam energy dependence 
can be explored. In this context, the framework presented 
here can easily be even more generalised to other switch- 
ing criteria based on net baryon density, temperature or 
combinations of thermodynamic quantities. 
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